// => extension of Matrix

// Includes

#include "matrix/matrix_ext.h"

#include <math.h>
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <sys/stat.h>

#include "matrix/matrix.h"

// Functions

/// @func: InitMatrix
/// >> init matrix with specific size and default value
/// @param: <row> the row of the matrix [ size_t ]
/// @param: <col> the column of the matrix [ size_t ]
/// @param: <val> the default value of the matrix [ double ]
/// @return: the matrix [ Matrix * ]
Matrix *InitMatrix(size_t row, size_t col, double val) {
  // get matrix pointer
  Matrix *mat = (Matrix *)calloc(1, sizeof(Matrix));
  // init
  mat->size[0] = row;
  mat->size[1] = col;
  mat->data = (double *)calloc(row * col, sizeof(double));
  // set all to value
  for (int i = 0; i < row * col; ++i) {
  }
  // return
  return mat;
}

/// @func: Ones
/// >> get a ones matrix with specific size
/// @param: <row> the row of the matrix [ size_t ]
/// @param: <col> the column of the matrix [ size_t ]
/// @return: a matrix fill with ones [ Matrix * ]
Matrix *Ones(size_t row, size_t col) { return InitMatrix(row, col, 1); }

/// @func: InitMatrixFromInput
/// >> init matrix from user input, default by row
/// @return: the matrix [ Matrix * ]
/// @descript: exit when the input is invalid
Matrix *InitMatrixFromInput(void) {
  // get row and column
  size_t row = 0, col = 0;
  printf("Please Enter the Size of the Matrix (row, col): ");
  scanf("%zu %zu", &row, &col);
  // get a raw matrix
  Matrix *mat = Zeros(row, col);
  // init matrix by input
  puts("Input Matrix Data:");
  for (size_t i = 1; i <= row; ++i) {
    for (size_t j = 1; j <= col; ++j) {
      double d = 0.0;
      scanf("%lf%*c", &d);
      SetMatrixValue(mat, i, j, d);
    }
  }
  // return matrix
  return mat;
}

/// @func: InitMatricesFromFile
/// >> init matrices from file
/// @param: <filename> the file path [ const char * ]
/// @param: <count> the MAXIMUM number of matrices [ size_t ]
/// @return: matrices [ Matrix ** ]
/// @descript:
///   * exit when file is not exist
///   * exit when content style is wrong
///   * exit when read file error
Matrix **InitMatricesFromFile(const char *filename, size_t count) {
  if (NULL == filename) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  }

  Matrix **mats = calloc(count, sizeof(Matrix *));
  size_t mcnt = 0;  // count the number of matrices
  // use read mode to open file
  FILE *fileHandler = fopen(filename, "r");
  // test whether file is opened successful
  if (!fileHandler) {
    perror("Open File");
    exit(EXIT_FAILURE);
  }
  // prepare buffer
  struct stat buf;
  if (stat(filename, &buf) == -1) {
    perror("Stat");
    exit(EXIT_FAILURE);
  }
  // get file content
  char *fileContent = malloc(buf.st_size);
  // read file
  while ((fscanf(fileHandler, "%[^\n] ", fileContent) != EOF) &&
         (mcnt < count)) {
    // find matrix head
    if (strcmp(">>", fileContent) == 0) {
      // get infomation about matrix
      if (fscanf(fileHandler, "%[^\n] ", fileContent) != EOF) {
        size_t row = 0;       // rows
        size_t col = 0;       // columns
        char orientation[4];  // orientation
        sscanf(fileContent, "%zu %zu %s", &row, &col, orientation);
        Orientation ort = ROW;  // default
        if (strcmp(orientation, "row") == 0) {
          ort = ROW;
        } else if (strcmp(orientation, "col") == 0) {
          ort = COL;
        }
        // get matrix data
        double *data = calloc(row * col, sizeof(double));
        size_t cnt = 0;
        char *elem = malloc(buf.st_size);
        while (fscanf(fileHandler, "%[^\n ] ", elem) != EOF &&
               cnt < row * col) {
          data[cnt++] = atof(elem);
        }
        // get matrix
        Matrix *mat = InitMatrixFromArray(row, col, data, ort);
        mats[mcnt++] = mat;
        free(elem);
        free(data);
      } else {
        perror("fscanf");
        exit(EXIT_FAILURE);
      }
    }
  }
  // clean
  free(fileContent);
  fclose(fileHandler);
  // adjust matrices size
  mats = realloc(mats, mcnt * sizeof(Matrix *));

  return mats;
}

/// @func: SaveMatrixToFile
/// >> save matrix data to file with text mode
/// @param: <filename> the file path [ const char * ]
/// @param: <mat> the matrix [ Matrix * ]
/// @descript: exit when `filename` or `mat` is null pointer
void SaveMatrixToFile(const char *filename, Matrix *mat) {
  if (NULL == filename || NULL == mat) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  }

  // open file
  FILE *fileHandler = fopen(filename, "a+");
  if (!fileHandler) {
    perror("Open File");
    exit(EXIT_FAILURE);
  }
  // start sign
  fputs(">>\n", fileHandler);
  // save matrix infomation, row by default
  fprintf(fileHandler, "%zu %zu row", mat->size[0], mat->size[1]);
  // save matrix data
  for (size_t i = 1; i <= mat->size[0]; ++i) {
    fputc('\n', fileHandler);
    for (size_t j = 1; j <= mat->size[1]; ++j) {
      fprintf(fileHandler, "%lf ", GetMatrixValue(mat, i, j));
    }
  }
  // end write
  fputs("\n<<\n", fileHandler);
  // close file
  fclose(fileHandler);
}

/// @func: PrintMatrixP
/// >> print matrix with a specific precision
/// @param: <mat> the matrix [ const Matrix * ]
/// @param: <prec> the precison [ int ]
/// @descript: exit when `mat` is a null pointer
void PrintMatrixP(const Matrix *mat, int prec) {
  if (NULL == mat) {
    fprintf(stderr, "mat is NULL!\n");
  }

  puts("matrix>>");
  for (size_t i = 1; i <= mat->size[0]; ++i) {
    printf(" [ ");
    for (size_t j = 1; j <= mat->size[1]; ++j) {
      printf("%*.*f ", prec + 4, prec, GetMatrixValue(mat, i, j));
    }
    printf("]\n");
  }
  puts("<<matrix");
}

/// @func: MatPow
/// >> the power of matrix
/// @param: <m> the matrix [ Matrix * ]
/// @param: <n> the power [ size_t ]
/// @return: the result [ Matrix * ]
/// @descript: exit when the matrix is not a square matrix
Matrix *MatPow(Matrix *m, size_t n) {
  if (m->size[0] != m->size[1]) {
    fprintf(stderr, "Not a Square Matrix Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *pM = Eye(m->size[0], m->size[1]);
  for (int i = 0; i < n; ++i) {
    Matrix *newPM = MatrixMul(m, pM);
    DestroyMatrix(pM);
    pM = newPM;
  }

  return pM;
}

/// @func: CatMatrix
/// >> concatenate two matrix vertically
/// @param: <m1> the first matrix [ const Matrix * ]
/// @param: <m2> the second matrix [ const Matrix * ]
/// @return: the result of concatenation [ Matrix * ]
/// @descript:
///   * exit when `m1` or `m2` is a null pointer
///   * exit when the size of `m1` and `m2` are not matched
Matrix *CatMatrix(const Matrix *m1, const Matrix *m2) {
  if (NULL == m1 || NULL == m2) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  } else if (m1->size[1] != m2->size[1]) {
    fprintf(stderr, "Unmatched Size Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *mat = Zeros(m1->size[0] + m2->size[0], m1->size[1]);

  for (int i = 1; i <= mat->size[0]; ++i) {
    for (int j = 1; j <= mat->size[1]; ++j) {
      if (i <= m1->size[0]) {
        SetMatrixValue(mat, i, j, GetMatrixValue(m1, i, j));
      } else {
        SetMatrixValue(mat, i, j, GetMatrixValue(m2, i - m1->size[0], j));
      }
    }
  }

  return mat;
}

/// @func: HorzCatMatrix
/// >> concatenate two matrix horizontally
/// @param: <m1> the first matrix [ const Matrix * ]
/// @param: <m2> the second matrix [ const Matrix * ]
/// @return: the result of concatenation [ Matrix * ]
/// @descript:
///   * exit when `m1` or `m2` is a null pointer
///   * exit when the size of `m1` and `m2` are not matched
Matrix *HorzCatMatrix(const Matrix *m1, const Matrix *m2) {
  if (NULL == m1 || NULL == m2) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  } else if (m1->size[0] != m2->size[0]) {
    fprintf(stderr, "Unmatched Size Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *mat = Zeros(m1->size[0], m1->size[1] + m2->size[1]);

  for (int i = 1; i <= mat->size[0]; ++i) {
    for (int j = 1; j <= mat->size[1]; ++j) {
      if (j <= m1->size[1]) {
        SetMatrixValue(mat, i, j, GetMatrixValue(m1, i, j));
      } else {
        SetMatrixValue(mat, i, j, GetMatrixValue(m2, i, j - m1->size[1]));
      }
    }
  }

  return mat;
}

/// @func: CopyMatrix
/// >> get a copy of matrix
/// @param: <mat> the matrix [ const Matrix * ]
/// @return: the copy of matrix [ Matrix * ]
/// @descript: exit when `mat` is null pointer
Matrix *CopyMatrix(const Matrix *mat) {
  Matrix *copy = Zeros(mat->size[0], mat->size[1]);

  for (int i = 0; i < mat->size[0] * mat->size[1]; ++i) {
    copy->data[i] = mat->data[i];
  }

  return copy;
}

/// @func: GetRowChangeMatrix
/// >> get a matrix that can change two rows
/// @param: <mat> the matrix [ const Matrix * ]
/// @param: <r1> the frist row [ size_t ]
/// @param: <r2> the second row [ size_t ]
/// @return: the row changed matrix [ Matrix * ]
/// @descript:
///   * exit when `mat` is a null pointer
///   * exit when row is out of boundary
Matrix *GetRowChangeMatrix(const Matrix *mat, size_t r1, size_t r2) {
  if (NULL == mat) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  } else if (r1 > mat->size[0] || r1 > mat->size[1] || r2 > mat->size[0] ||
             r2 > mat->size[1]) {
    fprintf(stderr, "Out of Boundary Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *rowChangeM = Eye(mat->size[0], mat->size[1]);

  SetMatrixValue(rowChangeM, r1, r1, 0);
  SetMatrixValue(rowChangeM, r1, r2, 1);
  SetMatrixValue(rowChangeM, r2, r2, 0);
  SetMatrixValue(rowChangeM, r2, r1, 1);

  return rowChangeM;
}

/// @func: InitLURes {PRIVATE}
/// >> get a lu decomposition result
/// @param: <lm> the matrix L [ Matrix * ]
/// @param: <um> the matrix U [ Matrix * ]
/// @param: <factor> the factor [ int ]
/// @return: the result [ LURes * ]
/// @descript: exit when `lm` or `um` is a null pointer
LURes *InitLURes(Matrix *lm, Matrix *um, int factor) {
  if (NULL == lm || NULL == um) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  }

  LURes *lures = calloc(1, sizeof(LURes));

  lures->factor = factor;
  // ! WARNING: may cause null pointer error
  lures->matrixL = lm;
  lures->matrixU = um;

  return lures;
}

/// @func: DestroyLURes
/// >> decomposition lu res
/// @param: <lures> the lu res [ LURes * ]
void DestroyLURes(LURes *lures) {
  if (NULL == lures) {
    return;
  }

  DestroyMatrix(lures->matrixL);
  DestroyMatrix(lures->matrixU);
  free(lures);

  lures = NULL;

  return;
}

/// @func: MatrixLUDecomposition
/// >> decomposite a matrix in LU form
/// @param: <mat> the matrix [ const Matrix * ]
/// @return: the result [ LURes * ]
/// @descript: exit when `mat` is a null pointer
LURes *MatrixLUDecomposition(const Matrix *mat) {
#define abs(x) ((x) < 0 ? -(x) : (x))
#define min(x, y) ((x) > (y) ? (y) : (x))
  // check prequest
  if (NULL == mat) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  }
  // init variables
  Matrix *copy = CopyMatrix(mat);
  Matrix *lMatrix = Eye(mat->size[0], mat->size[0]);
  int factor = 1;
  // the tolerance
  const double tol = 1e-5;
  // do decomposition
  for (size_t i = 1; i <= min(copy->size[0], copy->size[1]); ++i) {
    // check pivot whether a zero
    if (abs(GetMatrixValue(copy, i, i)) < tol) {
      // if pivot is a zero, then do row change
      for (size_t j = i + 1; j <= copy->size[0]; ++j) {
        if (abs(GetMatrixValue(copy, j, i)) > tol) {
          Matrix *rc = GetRowChangeMatrix(copy, i, j);
          // update Matrix L
          Matrix *newLM = MatrixMul(rc, lMatrix);
          DestroyMatrix(lMatrix);
          lMatrix = newLM;
          // update Matrix U
          Matrix *newM = MatrixMul(rc, copy);
          DestroyMatrix(copy);
          copy = newM;
          // update factor
          factor *= -1;
          // free row change matrix
          DestroyMatrix(rc);
          break;
        }

        if (j == copy->size[0]) {
          fprintf(stderr, "Bad Matrix Error!");
          exit(EXIT_FAILURE);
        }
      }
    }
    // do elimination
    for (size_t j = i + 1; j <= copy->size[0]; ++j) {
      if (abs(GetMatrixValue(copy, j, i)) > tol) {
        // get elimination matrix
        Matrix *elim = Eye(copy->size[0], copy->size[0]);
        SetMatrixValue(
            elim, j, i,
            -1 * GetMatrixValue(copy, j, i) / GetMatrixValue(copy, i, i));
        // update Matrix L
        Matrix *newLM = MatrixMul(elim, lMatrix);
        DestroyMatrix(lMatrix);
        lMatrix = newLM;
        // update Matrix U
        Matrix *newM = MatrixMul(elim, copy);
        DestroyMatrix(copy);
        copy = newM;
        // free elimination matrix
        DestroyMatrix(elim);
      }
    }
  }
#undef min
#undef abs

  return InitLURes(lMatrix, copy, factor);
}

/// @func: MatrixACofactor
/// >> get the algebric cofactor of matrix
/// @param: <mat> the matrix[ const Matrix * ]
/// @param: <row> the row [ size_t ]
/// @param: <col> the column [ size_t ]
/// @return: the algebric cofactor  [ double ]
/// @descript:
///   * exit when `mat` is a null pointer
///   * exit when `row` or `col` is out of boundary
double MatrixACofactor(const Matrix *mat, size_t row, size_t col) {
  if (NULL == mat) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  } else if (row < 1 || row > mat->size[0] || col < 1 || col > mat->size[1]) {
    fprintf(stderr, "Out of Boundary Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *subM = Zeros(mat->size[0] - 1, mat->size[1] - 1);
  size_t rc = 0;
  for (size_t i = 1; i <= subM->size[0]; ++i) {
    size_t cc = 0;
    if (i == row) {
      rc = 1;
    }
    for (size_t j = 1; j <= subM->size[1]; ++j) {
      if (j == col) {
        cc = 1;
      }
      SetMatrixValue(subM, i, j, GetMatrixValue(mat, i + rc, j + cc));
    }
  }

  return ((row + col) & 1 ? -1 : 1) * MatrixDeterminant(subM);
}

/// @func: InitAdjugateMatrix
/// >> get the adjugate of matrix
/// @param: <mat> the matrix [ const Matrix * ]
/// @return: the adjugate matrix [ Matrix * ]
/// @descript: exit when `mat` is a null pointer
Matrix *InitAdjugateMatrix(const Matrix *mat) {
  if (NULL == mat) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *adjM = Zeros(mat->size[1], mat->size[0]);

  for (size_t i = 1; i <= adjM->size[0]; ++i) {
    for (size_t j = 1; j <= adjM->size[1]; ++j) {
      SetMatrixValue(adjM, i, j, MatrixACofactor(mat, j, i));
    }
  }

  return adjM;
}

/// @func: LinearSolver
/// >> solve linear equations
/// @param: <A> the coefficient matrxi [ const Matrix * ]
/// @param: <b> the value vector [ const Matrix * ]
/// @return: the solution [ Matrix * ]
/// @descript:
///   * exit when `A` or `b` is null pointer
///   * exit when the size of `A` and `b` are not matched
///   * `b` must be column vector
Matrix *LinearSolver(const Matrix *A, const Matrix *b) {
  if (NULL == A || NULL == b) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  } else if (A->size[0] != b->size[0]) {
    fprintf(stderr, "Unmatched Size Error!");
    exit(EXIT_FAILURE);
  }

  return MatrixMul(InverseMatrix(A), b);
}

/// @func: JacobiLinearSolver
/// >> solve linear equations by Jacobi method
/// @param: <A> the coefficient matrix [ const Matrix * ]
/// @param: <b> the value vector [ const Matrix * ]
/// @param: <k> number of iterations [ size_t ]
/// @return: the solution [ Matrix * ]
/// @descript:
///   * exit when `A` or `b` is null pointer
///   * exit when the size of `A` and `b` are not matched
///   * `b` must be column vector
Matrix *JacobiLinearSolver(const Matrix *A, const Matrix *b, size_t k) {
  if (NULL == A || NULL == b) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  } else if (A->size[0] != b->size[0]) {
    fprintf(stderr, "Unmatched Size Error!");
    exit(EXIT_FAILURE);
  }

  Matrix *sol = Zeros(A->size[0], 1);
  // x_i^{k+1} = -\sum_{n=1, n != m}^r (D_{m m}^{-1} R_{m n} x_n^{k})
  //             + D_{m m}^{-1} b_m
  for (int i = 0; i < k; ++i) {
    for (int m = 1; m <= A->size[0]; ++m) {
      // the sum
      double sum = 0.0;
      for (int n = 1; n <= A->size[1]; ++n) {
        if (n != m) {
          sum += (1 / GetMatrixValue(A, m, m)) * GetMatrixValue(A, m, n) *
                 GetMatrixValue(sol, n, 1);
        }
      }
      SetMatrixValue(
          sol, m, 0,
          -sum + (1 / GetMatrixValue(A, m, m)) * GetMatrixValue(b, m, 1));
    }
  }

  return sol;
}

/// @func: GetEigSystem {# TODO #}
/// >> get matrix eigen value and eigen vectors
/// @param: <mat> the matrix [ const Matrix * ]
/// @return: the eigen system [ EigSys * ]
/// @descript: exit when `mat` is a null pointer
EigSys *GetEigSystem(const Matrix *mat) {
  if (NULL == mat) {
    fprintf(stderr, "Null Pointer Error!");
    exit(EXIT_FAILURE);
  }
  return NULL;
}
